Thermal metal in network models of a disordered two-dimensional superconductor 
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We study the universality class for localization which arises from models of non-interacting quasi- 
particles in disordered superconductors that have neither time-reversal nor spin-rotation symmetries. 
Two-dimensional systems in this category, which is known as class D, can display phases with three 
different types of quasiparticle dynamics: metallic, localized, or with a quantized (thermal) Hall con- 
ductance. Correspondingly, they can show a variety of derealization transitions. We illustrate this 
behavior by investigating numerically the phase diagrams of network models with the appropriate 
symmetry, and for the first time show the appearance of the metallic phase. 



PACS numbers: 73.20.Fz, 72.15.Rn 

The properties of quasiparticles in disordered super- 
conductors have been a subject of much recent interest. 
Within a mean field treatment of pairing, the quasiparti- 
cles are noninteracting fermions, governed by a quadratic 
Hamiltonian which may contain effects of disorder in 
both the normal part and the superconducting gap func- 
tion. Such Hamiltonians are representatives of a set of 
universality classes different from the three classes which 
are familiar both in normal disordered conductors and 
in the Wigner-Dyson random matrix ensembles. A list 
of additional random matrix ensembles, determined by 
these new symmetry classes, has been established rel- 
atively recently |jj. These additional random matrix 
ensembles describe zero-dimensional problems, and are 
appropriate to model a small grain of a superconduc- 
tor in the ergodic limit. In the corresponding higher- 
dimensional systems from the same symmetry classes, 
there can be transitions between metallic, localized, or 
quantized Hall phases for the quasiparticles [§-|J. The 
associated changes in quasiparticle dynamics must be 
probed by energy transport or (in singlet superconduc- 
tors) spin transport, rather than charge transport, since 
quasiparticle charge density is not conserved |2[. There 
are various possibilities for behavior, depending on the 
particular symmetry class considered. These have been 
studied theoretically using nonlinear sigma model meth- 
ods [gj, numerically J3|, and in quasi-one dimensional 
models || . An important question not addressed in such 
work so far, and which will not be considered here, is 
whether the self-consistent solution to the gap equation 
in the presence of disorder affects the universal statistical 
properties of the ensembles. 

In this paper we present extensive numerical results on 
a symmetry class with particularly rich phase diagram in 
two dimensions, class D. The symmetry may be realized 
in superconductors with broken time-reversal invariance, 
and either broken spin-rotation invariance (as in d-wave 
superconductors with spin-orbit scattering) or spinless or 



spin-polarized fermions (as in certain p-wave states) . The 
nonlinear sigma model for class D |lj] has been shown, in 
the two-dimensional case, to flow under the renormal- 
ization group to weaker values of the coupling constant 
|6|-|9|. The coupling constant is proportional to the in- 
verse of the thermal conductivity of the superconductor, 
and this flow implies that there is a phase in which there 
is a nonzero (indeed, diverging M) density of extended 
fcrmion eigenstates at zero excitation energy. A super- 
conductor described by this model would be in a thermal 
metal phase. We will refer to such a phase simply as a 
metallic phase. In addition, a phase with localized quasi- 
particles is a natural possibility, and — since time-reversal 
symmetry is broken — so is one with quantized Hall con- 
ductance. Our aim in the following is to investigate the 
appearance of these phases in simple models. 

As our starting point, we take versions of the net- 
work model pd| ] for a single-component fermion, which 
we specify in detail after first summarizing our findings. 
Disorder appears in the network model in the form of 
random scattering phases, and the symmetries of class D 
restrict scattering phases to the values and tt. Remark- 
ably, within this framework, different particular forms of 
disorder result in quite distinct physical behavior. We 
discuss three alternative choices. The first of these (CF) 
was introduced in work by Cho and Fisher |L2| with 
the intention of modeling the two-dimensional random 
bond Ising model (RBIM), which possesses a fermion 
representation with the symmetries of class D. In fact, 
as noted subsequently Jl^Jl^l j a precise mapping of the 
Ising model leads to a second version of the network 
model, which we label RBIM. In both these models, 
scattering phases with the value 7r appear in correlated 
pairs. A third model (also discussed in Ref. ||), which 
we denote by 0(1), arises naturally if one instead takes 
scattering phases to be independent random variables. 
Each model has two paranreters: a disorder concentra- 
tion, p (0 < p < 1), and a tunneling amplitude p|, a 
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FIG. 1. The phase diagram of the CF model obtained from 
our numerical calculations. 



(0 < a < 7r/2), which controls the value of the (thermal) 
Hall conductance at short distances. Our phase diagram 
for the CF model in the (p, a) plane is shown in Fig. fy. It 
contains a region of metallic phase, and two distinct local- 
ized phases, which can be identified with the ordered and 
disordered phases of the RBIM, or as regions with dif- 
ferent quantized Leduc-Righi (thermal Hall) conductivi- 
ties. As a consequence, three potentially different critical 
points occur: an insulator-to-insulator quantum-Hall- 
type transition; an insulator-to-metal transition; and a 
multicritical point at which all three phases meet. This 
phase diagram has the form proposed generically for class 
D in Ref. |j. In contrast, neither the RBIM nor the 
0(1) model supports all three phases: arguments that 
the metallic phase cannot appear in RBIMs with real 
Ising couplings are given in Ref. |L4|]; while in the 0(1) 
model we find no localized phase, in striking distinction 
to all network models studied previously. We show below 
how these differences can be understood by solving the 
models in one dimension and by considering them in two 
dimensions at weak tunneling. 

All these models represent coherent propagation of 
quantum-mechanical flux on a square lattice of directed 
links which meet at nodes, as illustrated in Fig. 0. Pla- 
quettes of the lattice can be divided into two sets, ac- 
cording to the direction of circulation around them. For 
general values of a, all plaquettes are coupled, but for 
a = the system separates into uncoupled plaquettes 
with clockwise circulation, while for a — n/2 it consists 
of uncoupled anticlockwise plaquettes. Disorder is in- 
troduced via a phase cj>i associated with each link, I. To 
make clear the constraints imposed in class D, recall that 
a Bogoliubov-de Gennes Hamiltonian with this symmetry 
may be written in terms of a purely imaginary Hcrmitian 
matrix jl| . The corresponding time evolution operator is 
purely real, restricting the generalized phase factors to 
be O(N) matrices for a model in which TV-component 
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FIG. 2. The network model. Values of the scattering ma- 
trix elements, ±cos(a) and ±sin(a), at nodes on each sub- 
lattice are indicated for p — schematically by ±C and ±S. 



fermions propagate on links, and to the values ±1 for 
N = 1, the case we treat. It is useful to consider the 
gauge-invariant total phase, modulo 2-7T, accumulated on 
passing around each elementary plaquette. In place of 
individual link phases, randomness can be characterised 
by the positions of flux lines which thread a subset of 
plaquettes, adding it to their phases. 

The models we study and the important distinctions 
between them are as follows. The CF model has transfer 
matrix tunneling parameters chosen negative with prob- 
ability p, and positive with probability 1 — p |12j . In con- 
sequence, flux lines appear in pairs at a node with proba- 
bility p: both members of a pair pass through plaquettes 
with the same circulation, but different pairs may belong 
to plaquettes with opposite circulation. The RBIM simi- 
larly has a pair of flux lines introduced at each node with 
probability p, but with the difference that all flux lines 
thread plaquettes of the same circulation |H],[l4| . Finally, 
the 0(1) model has link phase factors chosen negative 
with probability p and positive with probability 1 — p, 
so that the two members of a flux line pair are associ- 
ated with plaquettes of opposite circulation. Each of the 
models is invariant under the transformation p — * 1 — p, 
and so we consider only < p < 1/2. The CF and the 
0(1) models are both (statistically) self-dual for all p un- 
der a Kramers- Wannier transformation that takes a to 
7r/2 — a, leaving the line a = 7r/4 invariant. The RBIM is 
not self-dual, except at p — 0. Finally, the CF and 0(1) 
models are equivalent, under a gauge transformation, on 
the line p = 1/2. 

Some of the differences in the behavior of these 
three models can be illustrated by solving their one- 
dimensional versions, which consist of a single chain of 
links and nodes. In one dimension, disorder in the sign 
of the nearest-neighbor exchange interaction can be re- 
moved from the RBIM by gauge transformation, and the 



inverse localization length has the disorder-independent 
value £rb IM — arc tanh(| sin(a)|), finite for all a^O, n/2. 
An elementary calculation gives for the CF model 



£CF = |1-2P|£ 



RBIM ' 



(1) 



so that £cf diverges as p — * 1/2 but is otherwise finite, 
while for the 0(1) model Com = for all p ^ 0,1. 
Thus, the localization properties of the one-dimensional 
CF model at p ^ 1/2 are like those of models belonging 
to the Wigner-Dyson universality classes, in that states 
are localized, while the absence of localization in the 0(1) 
model mirrors that found previously in quasi-one dimen- 
sional class D systems M. 

A second useful approach illustrating differences be- 
tween these models is to consider their two-dimensional 
versions in the limit of weak inter-plaqucttc tunneling 
{a <Sil or 7r/2 -a<l) and weak disorder (p <C 1). We 
do this in terms of the discrete-time evolution operator 
U in a closed system, which is an Ni x Ni unitary matrix 
for a network of Ni links [ p~5|Jl6| ; a similar approach to 
a different free fermion representation of the RBIM was 
used in Ref . |17]] . The eigenvalues of U lie on the unit cir- 
cle and may be written as e~ le . The real e (— tt < e < it) 
play the role of excitation energy eigenvalues, and are 
distributed symmetrically in pairs around e = because 
U is a real orthogonal matrix. Long-time properties are 
determined by the part of the spectrum near e = 0, on 
which we now focus. At zero tunneling, it is sufficient 
to examine an isolated plaquette. In our disorder-free 
reference system, the evolution operator for a single pla- 
quette satisfies U 4 = — 1, and hence e = ±7r/4, ±37r/4. 
For a single plaquette with a flux line added, U 4 = 1 
and e = 0, n, ±7r/2. Turning on weak tunneling, it is 
clear that the spectrum near e = for a large system will 
arise by hybridisation of the e = states from plaque- 
ttes with flux lines. In both the RBIM and CF models, 
there are two scales for this hybridisation, because flux 
lines appear in the system in adjacent pairs associated 
with plaquettes of the same circulation. The first conse- 
quence of tunneling is to remove the degeneracy within 
each pair, yielding approximate eigenvalues e = ±eo- At 
small p, pairs are dilute and tunneling between differ- 
ent pairs is not sufficient to generate extended states at 
e = 0. By contrast, for the 0(1) model in this regime, 
there is only one scale for hybridisation, since single flux 
lines appear independently on the set of weakly-coupled 
plaquettes. As a result, metallic behavior is not excluded 
even at p, a <C 1 . 

Our results from numerical simulation supplement this 
qualitative discussion. We study the CF and 0(1) models 
in cylindrical geometry via the transfer matrix T, obtain- 
ing the positive Lyapunov exponents, < i>\ < . . . < vm 
in a system of width M' = 2M links. A crucial tech- 
nical aspect of these calculations is our discovery that 
the standard algorithm |l9LpO| has a serious instability to 



roundoff errors throughout much of the phase diagram of 
both models. More specifically, we find that the smallest 
positive Lyapunov exponent, v\, may be either identi- 
cally zero or exceptionally small {i>\ <§; M _1 ). (The first 
happens in the 0(1) model for all p and a, and in the 
CF model on the self-dual line a — 7r/4; the second hap- 
pens in the metallic phase of the CF model.) Under these 
circumstances, numerical noise from roundoff errors gen- 
erates a systematic positive error in the value obtain for 
v\ . From an analytical theory p8[ of the instability, we 
find that the error in v\ decreases with reduced noise 
amplitude r\ only as |log(7y)| _1 . This instability can be 
cured by making explicit use in numerical calculations of 
the structure imposed on T by current conservation and 
the symmetry of class D. 

In detail, T has the polar decomposition 



T = 



A t 
A 2 



cosh(7) sinh(7) 
sinh(7) cosh(7) 



A? 





A\ 



(2) 



where A\ . . . A& are M x M real orthogonal matrices and 
7 is an M x M real diagonal matrix. It follows that T T T 
is diagonalized by the transformation B T T T TB, where 



B 



A 3 A 3 
Ai -At 



(3) 



The standard method for calculating Lyapunov expo- 
nents numerically involves acting with the transfer ma- 
trices for successive slices of the system on a set of M or- 
thogonal vectors, and reimposing orthogonality by means 
of Gram-Schmidt transformations JT^^y] . If all Lya- 
punov exponents are separated by gaps, this set of vectors 
converges to the eigenvectors of T T T associated with the 
first M exponents. Convergence rates are determined by 
the sizes of gaps between successive exponents. In the 
present case, convergence rates are seriously reduced if 
v\ approaches zero, so that the gap between the smallest 
positive and largest negative exponents vanishes. More- 
over, numerical noise ultimately limits the extent of con- 
vergence, and leads to an erroneously large value for v\. 
To overcome this, we impose on the M vectors concerned 
not simply orthogonality but instead the fact that their 
first M components separately form an orthogonal ma- 
trix A 3 , and their last M components form A±, as is 
evident from Eq.|3j. The results we obtain in this way for 
the CF model differ significantly from those of Ref. J 12| . 
Evidence in support of the phase diagram of Fig.jl for 
the CF model is presented in Fig. [| On the self-dual line 
(a = 7r/4) we believe that v\ is identically zero (as in the 
one-dimensional 0(1) model). For example, at p = 1/2 
and a = tt/A we obtain in systems of length L = 5 TO 5 the 
bounds v x < 1.5-1CT 3 at width M' = 4 and v x < 1.5-1CT 4 
at M' — 256. In order to search for a possible multi- 
critical point on the self-dual line, we therefore examine 
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FIG. 3. Behavior of the CF model in systems of width 
W = 64 (o), 128 (□), and 256 (O). (a) The self-dual line: 
M'f2 as a function of p. (b) Quantum-Hall-type transition: 
M'v\ as a function of sin 2 (a) at p = 0.1. (c) Insulator - metal 
transition: M'v\ as a function of p at sin 2 (a) = 0.19. 



the behavior of V2 El] . If there is a multicritical point 
at p = pmCi one expects the amplitude ratio M'v-x to 
show three regimes at large M', as a function of p. For 
V < Pmc, scaling flow is towards smaller p and M'i>i has 
a p-independent value governed by the critical point at 
p = 0. At p — PmCi a distinct limiting value arises from 
the multicritical point. And for p > j>mc, scaling flow 
of the conductivity in the metallic phase towards larger 
values means that M'v% will slowly decrease towards zero 
with increasing M'. The data shown in Fig.|3|a are con- 
sistent with this scenario, although the position of the 
multicritical point is not well-determined: we find the 
bounds 0.05 < pmc < 0.15. A quantum-Hall-type tran- 
sition is observed on crossing the self-dual line by varying 
a at hxed p < pmc, a s illustrated in Fig.gb: M'v\ in- 
creases with M' for a 7^ 7r/4 (localization) and vanishes 
as a — > 7r/4 (delocalization). (This transition is expected 
J7| to be in the universality class of the pure Ising tran- 
sition, because the disorder strength scales towards zero, 
as in the RBIM at small p.) We determine the position of 
the metal-insulator phase boundary from the variation of 
M'v\ with p and M' at fixed a ^ 7r/4, shown in Fig.||c. 
Since M'v\ decreases rapidly with increasing M' in the 
metal and increases with M ' in the insulator, the critical 
point, pc(oi), is identified by the crossing of curves for 
different M' . In this way, we arrive at the phase diagram 
for the CF model displayed in Fig. [I]. 

We believe that the 0(1) model has only a metallic 
phase, and has v\ identically zero for all p ^ 0. Our 
calculations cover the range 0.1 < p < 0.5 and 0.1 < 
sin 2 (a) < 0.5. If the model were to support a localized 
phase, it should appear at small p, a. As an illustration 
of the absence of localization, at p — 0.1, sin 2 (a) =0.1, 



we calculate for M' = 16: v\ < 10~ 3 in the 0(1) model, 
while v\ = 0.83 in the CF model. 

In summary, we find that two-dimensional models for 
localization in the symmetry class D can have quite dif- 
ferent behavior according to the form of disorder. Several 
additional points deserve emphasis. The metallic phase 
of the CF model is self-dual, as is its multicritical point. 
By contrast, the RBIM is not self-dual but has higher su- 
persymmetry at its multicritical point jl(| . There is little 
reason to suppose that these two multicritical points are 
in the same universality class. Separately, the apparent 
absence of an insulating phase in the 0(1) model is re- 
markable, because the bare conductivity becomes small 
when a — > or n/2. Recently, it has been emphasized 
that the target manifold of the class D nonlinear sigma 
model is not connected p], and this means that domain 
wall excitations can occur in the sigma model, which 
must be described by additional parameters, and have 
not been taken into account in weak-coupling analyses 
so far. It is likely that these domain walls in the sigma 
model language are connected with the richness of phases 
in this symmetry class. In that context, the 0(1) model 
with p = 1/2 is known to be a special case, since it maps 
to a sigma model without domain walls: this fact suggests 
that proliferation of domain walls may be necessary for 
localization §0Q. 
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